/* ------------------------------- Gas-phase ------------------------------- */
Foam::Info << "\nInitialize field variables of the gas-phase" << Foam::endl;

// Density
rho = p / Rg / T;
rho.correctBoundaryConditions();

// Properties
for (int i = 0; i < mesh_gas.nCells(); i++)
{
    calc_gas_property(T[i], lambda[i], mu[i], Cp[i]);
}
mu.correctBoundaryConditions();
Cp.correctBoundaryConditions();
lambda.correctBoundaryConditions();

// Signed-distance
for (int i = 0; i < pointMesh_gas.size(); i++)
    phi_gas[i] = mesh_gas.points()[i].z() - plane_z0;
phi_gas.correctBoundaryConditions();

for (int i = 0; i < mesh_gas.nCells(); i++)
{
    phi_C_gas[i] = mesh_gas.C()[i].z() - plane_z0;
    calcGrad(mesh_gas, phi_gas, i, grad_phi_gas[i]);
}
phi_C_gas.correctBoundaryConditions();
grad_phi_gas.correctBoundaryConditions();

// Immersed-Boundary flags
identifyIBCell(mesh_gas, phi_gas, phi_C_gas, cIbMarker);
updateIBMask(cIbMarker, cIbMask);
cIbMarker.correctBoundaryConditions();
cIbMask.correctBoundaryConditions();

// Volume flux (createPhi.H)
USn = Foam::fvc::interpolate(U) & mesh_gas.Sf();
setBdryVal(mesh_gas, U, USn);

// Mass flux (compressibleCreatePhi.H)
rhoUSn = Foam::fvc::interpolate(rho*U) & mesh_gas.Sf();
setBdryVal(mesh_gas, rho, U, rhoUSn);

/* ------------------------------ Solid-phase ------------------------------ */
Foam::Info << "\nInitialize field variables of the solid-phase" << Foam::endl;

// Density, Temperature & Properties
for (int i = 0; i < mesh_solid.nCells(); i++)
{
    T_solid[i] = T0;
    if (mesh_solid.C()[i].x() < 0.3*L || mesh_solid.C()[i].x() > 0.7*L)
    {
        rho_solid[i] = rho_AP;
        c[i] = c_AP;
        lambda_solid[i] = lambda_AP;
    }
    else
    {
        rho_solid[i] = rho_HTPB;
        c[i] = c_HTPB;
        lambda_solid[i] = lambda_HTPB;
    }
}
T_solid.correctBoundaryConditions();
rho_solid.correctBoundaryConditions();
c.correctBoundaryConditions();
lambda_solid.correctBoundaryConditions(); 

// Signed-distance
for (int i = 0; i < pointMesh_solid.size(); i++)
    phi[i] = mesh_solid.points()[i].z() - plane_z0;
phi.correctBoundaryConditions();

for (int i = 0; i < mesh_solid.nCells(); i++)
{
    phi_C[i] = mesh_solid.C()[i].z() - plane_z0;
    calcGrad(mesh_solid, phi, i, grad_phi[i]);
}
phi_C.correctBoundaryConditions();
grad_phi.correctBoundaryConditions();
